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We present a new likelihood-ratio ranking statistic for use in searches for gravitational waves 
from the inspiral and merger of compact object binaries. Expanding on previous work, the ranking 
statistic incorporates a model for the correlations in the signal-to-noise ratios with which signals will 
be seen in a network of ground-based antennas while retaining an algebraic procedure for mapping 
ranking statistic values to false-alarm probability. Additionally, the ranking statistic enables the im¬ 
plementation of a rigorous signal rate estimation technique. We implement the ranking statistic and 
demonstrate its use including signal rate estimation in an analysis of a simulated signal population 
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I. INTRODUCTION 

An essential ingredient of searches for gravitational 
waves is a statistic by which to rank candidates from least 
signal-like to most signal-like. The likelihood ratio (LR) 
has been shown to provide the most powerful detection 
statistic at fixed false-alarm probability [1], but due to 
technical challenges all ranking statistics used to date in 
searches for compact binary coalescences (CBCs) have 
been approximations of a LR. To be useful, the rank¬ 
ing statistic must not only be effective at ranking signals 
above noise, but it must not be computationally costly 
to implement, and there must be a known mapping from 
ranking statistic value to false-alarm probability. 

In [2], Biswas et al. developed theoretical aspects of 
the use of the LR [3, Section 7.6] as the ranking statistic 
for searches for gravitational-waves (GWs) from CBCs. 
For example, they showed that when a LR-based rank¬ 
ing statistic is used the detection efficiency is not dimin¬ 
ished when the volume of the parameter space over which 
a search is conducted is increased, in contrast to tech¬ 
niques using other ranking statistics [4]. That work was 
preceded by an earlier effort to incorporate the Virgo 
antenna [5] into a search for CBCs [6], wherein the sub¬ 
stantial difference in sensitivity of the Virgo antenna next 
to the three LICO antennas [7] of the day motivated the 
need to explicitly account for the relative sensitivities of 
all antennas in the ranking statistic. Biswas et al. demon¬ 
strated an implementation of a LR ranking statistic on 
LICO’s S4 data, but it was implemented as a rescaling 
of the traditional “combined effective signal-to-noise ra¬ 
tio (SNR)” ranking statistic [2, equation (26)], and so 
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while the LR allowed events from times when different 
instruments were operating or from different regions of 
mass parameter space to be ranked on a common scale it 
otherwise provided no information not already contained 
in the combined effective SNR, and they continued to 
rely on time-shift analyses [8] to map LR values to false 
alarm probability (FAR). In the related work [9], Biswas 
et al. demonstrate the use of the LR to combine the re¬ 
sults from several, possibly degenerate or uninformative, 
searches for GWs into a single unihed search result. 

In [10], Aasi et al. built upon the histogram-based LR 
ranking statistic described in [11] for coincidences in a 
time-frequency burst search to develop a LR-like ranking 
statistic for searches for GW bursts from cosmic string 
cusps. Subsequently, in [12] Cannon et al. expanded upon 
that technique to develop an histogram-based LR-like 
ranking statistic for searches for GWs from CBCs, and 
demonstrated its use. There has also been work to ap¬ 
proximate LR ranking statistics using space-partitioning 
decision trees in searches for GW bursts [13], and in CBC 
searches [14]; in the latter the mapping from ranking 
statistic to FAR was accomplished through the use of 
time-slides while in the former the ranking statistic was 
applied to gamma-ray burst (GRB)-triggered searches 
and the mapping was constructed using on-source/off- 
source comparisons. In contrast, the ranking statistic 
developed in [12] had the property that if applied to 
the candidates collected from a search pipeline possess¬ 
ing certain simplifying characteristics, it was possible to 
compute the mapping from ranking statistic to FAR al¬ 
gebraically. 

The ranking statistic in [12] was an approximation of 
the LR in which it was assumed that the characteristics 
of a signal-like event in one instrument were indepen¬ 
dent of that same event’s characteristics in other instru¬ 
ments — that the probability density function (RDF) in 
the ranking statistic’s numerator could be factored into a 
product of single-instrument terms. This approximation 
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simplified the software needed to implement the rank¬ 
ing statistic and afforded a simple numerical integration 
scheme for computing the FAP for any value thereof. 

Although the ability to compute FAPs algebraically 
was a significant improvement over previous techniques, 
three limitations remained, (i) The ranking statistic as¬ 
sumed the SNRs at which a signal is observed in the var¬ 
ious instruments are independent random variables. For 
low SNRs where plausible early detections are expected 
this is a good approximation but for higher SNRs this 
approximation yields a less-than-optimal ranking statis¬ 
tic in that it does not properly penalize candidates with 
implausible SNR combinations that have resulted from, 
e.^., non-stationary noise in one detector, (ii) No at¬ 
tempt was made to include within the ranking statistic 
knowledge of which instrument combination had identi¬ 
fied the candidate, leading to a unique ranking statistic 
scale for each combination of instruments, something the 
work of Biswas et al. had already addressed, (iii) Finally, 
although the FAP for each value of the ranking statistic 
could be computed, with only an approximate PDF for 
the numerator a true-alarm probability (TAP) (proba¬ 
bility of at least one candidate at or above a given value 
of the ranking statistic given the presence of a signal) 
could not be computed which excluded the possibility of 
combining the ranking statistic with the rate estimation 
technique of Farr et al. [15]. 

Here, we address these limitations by showing how to 
compute and incorporate the probability of observing sig¬ 
nals and noise events in various combinations of instru¬ 
ments into both the numerator and the denominator, and 
how to incorporate the joint PDF for the SNRs in the nu¬ 
merator. These modifications will require us to abandon 
the use of explicit numerical integration for computing 
the ranking statistic PDF in signal-free data — the first 
step in mapping ranking statistic values to FAPs — and 
switch to an importance-weighted sampling procedure, 
the details of which will also be explained. 


II. EXAMPLE 

To demonstrate the application of the ranking statistic, 
we have performed a test analysis of about 5.23 x 10® s 
of stationary Gaussian noise simulating the Advanced 
LIGO and Virgo antennas and possessing the spectral 
densities described as the “early” noise curves in [16, 17]. 
A population of binary neutron star (BNS) merger events 
was injected into the data from a population of sources 


distributed uniformly in volume to a distance of 300 Mpc, 
uniformly in component masses in the range 1 M0-3 Mq, 
with Gaussian-distributed dimensionless spins centred on 
0 with a standard deviation of 0.4 a cut-off of 0.7 and 
isotropic orientations, and occurring at a mean rate of 
10 Mpc“®Ma“^. This is not meant to be astrophysically 
realistic [18], but to provide a useful sample of detectable 
signals. Gomparing the injected waveforms to the noise 
spectra and accounting for the relative orientations of 
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FIG. 1. Expectation value for the number of coincident candi¬ 
dates to be recovered from the signal population injected into 
test data set as a function of the single-detector SNR thresh¬ 
old applied in the search assuming an SNR recovery efficiency 
of 0.975. 

the antennas and the sources at the time of each sim¬ 
ulated event, one can determine the SNR at which the 
signal should be seen in each antenna, and from that es¬ 
timate how many signals should be seen in two or more 
antennas (z.e., form a coincident candidate event) as a 
function of the single-detector SNR threshold applied. 
This is shown in FIG. 1. The data was filtered using 
a template bank consisting of 6278 non-spinning BNS 
waveforms with component masses in the range 1 Mq- 
3Mq, and having a minimal match of 0.975 (with re¬ 
spect to the template family). The template bank was 
laid out using the lalapps_tmpltbank program from the 
LALSuite software package [19], and candidate events 
collected using the gstlal_inspiral program from the 
GstLAL software package [20, 21]. The plots and results 
shown in what follows are taken from this analysis. 

III. RANKING STATISTIC, C 

The ranking statistic is the LR 


^ ({-Dhhi, ^hli, • ■ •} > {HI, LI,...} , Phi, Xhi7 dLi, Xm> ■ • ■ , 

= C( IglTfg) -P({-Phhi,-Dhli,---},{H1, LI,...}, PHI, xgii,pLi,xL^---l^^ signal) ^ 
H ({Dhhi,-Dhli, • • , (HI, LI,...} , Phi, Xhi’Pli, xL> ■ ■ ■ noise) 
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where {-Dhhi? -Dhlij ■ • ■} is the set of horizon distances 
for all instruments in the network at the time the event 
is observed (see Appendix A), {HI,LI,...} is the spe¬ 
cific set of instruments in which the event has been ob¬ 
served, p and are the template matched-filter SNRs 
and values [22, and references therein] for the can¬ 
didate in each of those instruments, and 9 are the in¬ 
trinsic parameters of the template. The time at which 
the event is observed enters implicitly through the hori¬ 
zon distances which fluctuate in time as environmen¬ 
tal factors and operating conditions at the observato¬ 
ries change; in fact, {L>hhi 7 ^hlI) ■ • ■} is best thought 
of as simply a clock, and since time influences the dis¬ 
tributions of other parameters through their dependence 
on the (time-dependent) instrument sensitivities, the la¬ 
belling of times on the clock is most conveniently done 
with the instrument sensitivities themselves. We’ll say a 
few more words about the time dependence and the rank¬ 


ing statistic below. The template’s intrinsic parameters, 
9, can be a label identifying a template in the search, or 
identifying a group of templates if that is found to be 
convenient. Anticipating what follows, we have factored 
9 out of the PDFs on the right-hand side of (1). For 
brevity we will drop 9 from the notation in the remain¬ 
der of Section III; it should be remembered that we will 
be deriving PDFs for a choice of 9. 

The differences between the times at which the event 
is observed at each location on Earth are not used, nor 
are the phases with which the waveforms are recovered 
at each observatory. These omissions are an obvious av¬ 
enue for future improvements but, for signals, the inter¬ 
antenna time delays and waveform phases are correlated 
with the SNRs and their inclusion would make the nu¬ 
merator substantially more difficult to compute numeri¬ 
cally. 

We factor the numerator of the fraction in (1) as 


= P (I-Dhhi, ^HLi7 ■ ■■})P ({HI, LI,...}[ {Dhhij-Dhli, ■ • ■} , signal) 

X R(pHi,PLi,---|{A'HHi7^HLi,---},{Hl,Ll,...},signal) P (xLt I Anst, signal) . (2) 


That is, we write it as the probability of the collection of 
instruments having some given sensitivities multiplied by 
the probability that a (detectable) signal is visible in a set 
of those instruments (and no others) given the sensitivi¬ 
ties of all instruments multiplied by the joint probability 
density of the SNRs in that set of instruments multi¬ 
plied by the product of the probability densities for the 
X^ values given the SNRs and the instruments. In par¬ 
ticular, we assume here that except for their correlation 
with SNR the x^ values observed in different instruments 
for the same event are statistically independent random 
variables, i.e., that the residual from which the x^ is 
computed is dominated by instrumental noise and not 
by systematic waveform errors (which would be common 
to all instruments in the network). The validity of this 
assumption can be checked visually with scatter plots like 
those shown in FIG. 2. 

We assume the noise processes in the different instru¬ 
ments are statistically independent of each other and in¬ 
dependent of the sensitivity of the antennas to GWs and 
so the denominator of the fraction in (1) can be factored 
into several low-dimensional terms. 

= P ({Dhhi, ■■■})P ({HI, LI,...} [noise) 

X Yl P(pi„st,xLt|noise) . (3) 

instG{Hl,Ll,...} 

The factor P ({Dhhi, ^hlI) ■ ■ •}) is common to the 
numerator and denominator and, so, for the purpose of 
evaluating C the horizon distances appear only paramet¬ 


rically and only in the numerator. Later, when using the 
numerator and denominator to map C values to FAP and 
TAP the factor P ({Phhi, -DhlI) ■ • •}) will be included in 
those integrals. Since we understand {Phhi: ^hlii • ■ •} 
to be the labelling of a clock, P ({Phhi,-Dhli) • ■ •}) = 
I/livetime. In the following subsections we show how 
each of the remaining terms in (2) and (3) is obtained. 

A. P ({HI, LI,...} I {Dhhi, Phli, . • •} , signal) 

To obtain the probability that a signal yields coinci¬ 
dent above-threshold events in exactly some combination 
of instruments, we begin by assuming that the coinci¬ 
dence time window is sufficiently large that if a signal 
is seen above the SNR threshold in two instruments the 
threshold crossings will occur within that pair’s coinci¬ 
dence window with certainty. This assumption is essen¬ 
tially equivalent to requiring the coincidence window for 
each pair of instruments to be at least as large as the 
sum of their reciprocal bandwidths plus the light travel 
time between them (at most few tens of milliseconds for 
ground-based laser interferometer antennas). If a more 
sophisticated ranking statistic is considered in the future 
in which the time-of-arrival of the events in each of the 
instruments is treated more carefully than to simply con¬ 
struct a pass/fail test, or if the coincidence window for 
the pass/fail test is small, then this assumption would 
need to be revisited. 

Given our assumption about the formation of coinci- 
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FIG. 2. Statistical independence of x^ observed in different 
instruments. Plots show reduced x^ values from distinct (but 
irrelevant) instruments in coincidences arising from software 
injections subject to the constraint that the observed SNRs 
be in the given ranges. There is a trend towards larger x^ 
values as the SNR of the injection is increased, but at a given 
SNR the values in different instruments are not signifi¬ 
cantly correlated indicating that instrument noise dominates 
the residuals over template/signal mismatch. Note that this 
is despite the search being conducted with non-spinning tem¬ 
plates while the injections have spin and include orbital pre¬ 
cession. 


dences, we need only consider the probability that the 
signal is seen above the SNR threshold. The probabil¬ 
ity that a detectable signal is visible in exactly some 
collection of instruments, {HI,LI,...}, (and no oth¬ 
ers) can be obtained to satisfactory accuracy by Monte 
Carlo integration. Referring to Appendix A, a uniformly- 
distributed source position is selected on the sky and the 
antenna responses to the -I- and x polarizations com¬ 
puted, the cosine of the orbit inclination is drawn uni¬ 
formly from [—1,-|-1], and from the antenna responses 
and orbit inclination the quantity D is computed for 
each instrument in the network using (A4). The num¬ 
ber of sources visible above the SNR threshold p* in two 
instruments whose sensitivities are characterized by Di 
and D 2 is 


oc 





(4) 


while the number of sources visible to two instruments 
and not a third is 






dD, 


(5) 

and so on. All of these integrals diverge for the same 
reason that (A6) diverges: Qi yf 0 as D —)■ 00 . The 
probabilities that we seek will be the ratios of sums of 
integrals like (4) and (5), and even though the integrals 
diverge it might be the case that the ratios we seek are 
well-defined; or it might be that the ratios can be made to 
be well-defined by introducing a regularization, possibly 


one derived from the constraint that we seek a ranking 
statistic to extremize detection efficiency at fixed false 
alarm rate (FAR). 

Unfortunately we have not been able to make progress 
in this direction. Instead, we proceed as follows. We 
take (A7) to give the rate (up to an irrelevant constant) 
at which signals in the given direction and with the given 
polarization and given orbit inclination are seen by each 
instrument in the network, and then assume that all sig¬ 
nals visible to the least sensitive instrument are visible 
to the next least sensitive instrument with certainty, and 
so on. This works best if the instruments are not all very 
similar in sensitivity, which is likely to be the case as the 
network of ground-based GW antennas is commissioned 
in a staggered schedule. With that assumption, from 
differences and ratios of the rates one can compute the 
probability that a source will be visible to some combi¬ 
nation of instruments (and no others) assuming that it is 
visible to at least two. These probabilities are added to a 
histogram and the process of selecting a sky location and 
so on is repeated. When sufficient iterations have been 
made, the histogram is normalized. We find that after 
500,000 iterations the variance of the result is reduced to 
one part in 0(10®), and a Python implementation of the 
sampling loop using antenna response code implemented 
in C in LALSuite [19] completes in about 30 s on typical 
modern hardware. An example for p* = 4 is shown in 
the bottom-left panel of FIG. 3. 


B. 

P(pHi,PLi,... I {Uhhi.T'hli, • • •} , {HI, LI,...} , signal) 

The calculation of the joint PDF of the SNRs involves 
the same framework as was used in Section III A, and 
encounters some of the same difficulties. The PDFs are 
computed using Monte Carlo integration of the density 
at each site in a multi-dimensional SNR histogram. Each 
axis is binned using an tan“^ In binning as described in 
Appendix B. Referring to Appendix A, a uniformly- 
distributed direction is selected on the sky and the an¬ 
tenna responses to the -I- and x polarizations computed, 
the cosine of an orbit inclination is drawn uniformly from 
[—1, -1-1], and from these the quantity D is computed for 
each instrument in the network using (A4). 

We identify the most sensitive of the instruments that 
must see the source, and step through a sequence of bins 
of nominal SNR, pQ, in that instrument. Each bin has 
a lower bound, an upper bound, and a central value. 
Erom the D’s and each bin’s central value, the po in 
all other instruments are computed. With respect to 
the instruments that must not see the signal, we make 
the assumption that if the nominal SNR in those instru¬ 
ments is below the detection threshold, p*, then they are 
all blind to it with certainty, and if above p* in one of 
them it sees it with certainty. We are approximating the 
Marcum Q-function with a Heaviside function for per¬ 
formance purposes, and will rely on the superposition of 
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computing a vector of Rician-distributed random vari¬ 
ables (RVs) having noncentrality parameters given by 
the bin centre’s co-ordinates to give a vector of observed 
SNRs, and adding — Df to that location in the binned 
PDF. A new sky location, etc., is chosen and the whole 
process is repeated. Note that this procedure results in 
each bin in the SNR PDF containing a number propor¬ 
tional to the rate of signals in that bin, not the PDF in 
that bin. 

After a satisfactory number of iterations has been com¬ 
pleted, the array of signal rates is convolved with a Gaus¬ 
sian density estimation kernel, all bins below p* in any 
instrument are zeroed, and the array divided by its sum 
(so that its sum is 1). Finally each bin is divided by its 
volume in SNR to yield the properly-normalized PDF. 
Examples of PDFs obtained this way are shown in FIG. 
5. For these, the tan“^ln binnings along each axis are 
given by x\o = 3.6, a;hi = 120, n = 100, 80,000 itera¬ 
tions of the outer sky-location loop were performed, and 
a Gaussian density estimation kernel with a standard de¬ 
viation of 1.875 bins was employed. 


FIG. 3. (Top Left) Total number (in parentheses) of non¬ 
coincident, “noise”, events observed in each instrument and 
percent of total. (Top Right) Total number (in parenthe¬ 
ses) of noise coincidence events predicted to be observed for 
each combination of participating instruments, and the per¬ 
cent of total. This is P ({HI, LI,...} |noise), see Section III D. 
(Bottom Right) Total number (in parentheses) of coincidence 
events observed for each combination of participating instru¬ 
ments, and the percentage of the total. Notice the excel¬ 
lent agreement with the predicted percentages. (Bottom Left) 
Probability that a signal visible to at least two instruments 
is visible to the given combinations of instruments. This is 
P ({HI, LI,...} I {Dhhi, Dhli, • • •} , signal), see Section HI A. 



FIG. 4. The number of sources in the plug of space is propor¬ 
tional to its volume, which is oc Df — Df. 


many trials combined with a convolution of the histogram 
with a density estimation kernel to hide the approxima¬ 
tion. The binning of po in tlie most sensitive instrument 
is started at 1, and terminated when this latter condition 
is met. Using (A4) and D in the most sensitive instru¬ 
ment, the lower and upper bounds of each po bin can be 
converted into upper and lower bounds, respectively, on 
the physical distances of the corresponding sources. The 
number of sources in the bin is proportional to the dif¬ 
ference of the cubes of these, see FIG. 4. The SNR PDF 
is computed by iterating over the po binning, at each bin 


C. P (x^I p, signal) 

We begin by obtaining P (p, x^/p^| signal) as described 
in [12], using an analytic expression obtained by assum¬ 
ing recoverable signals are to be found in glitch-free, 
Gaussian, noise. As before, the implementation samples 
the PDF on a discrete rectangular grid of x^/p^ vs. p in 
order to better fit the natural shapes of the PDF’s con¬ 
tours to the rectangular grid of bins. For both the p and 
X^/p^ axes we now use the tan“^ In binnings described 
in Appendix B. The p axes are binned using x\o = 3.6, 
cchi = 70, n = 260; the jp^ axes using x\o = .001, 
Xhi = .5, n = 200. 

Once that PDF is obtained, columns of constant p are 
renormalized so that their integrals are 1, converting the 
function represented by the array from a two-dimensional 
density to a 1-dimensional density in x^/p^ parameter¬ 
ized by p. 

D. P({Hl,Ll,...}|noise) 

The search algorithm employs a bank of two-phase 
matched filters applied to the strain time series recorded 
from each of the instruments. When the magnitude of 
the output of a filter crosses a preset threshold, p*, a 
peak finder identifies the sample with the highest SNR, 
and that point in that filter’s output is recorded as an 
event, (sometimes called a “trigger”, borrowing language 
from particle experiments). A “coincidence” occurs when 
two or more instruments register events from the same 
filter within some window of time, \ti — t 2 \ < ti 2 = 
5t + \xi — X 2 \ /c, where xi and X 2 are the positions of 
the two antennas and c is the speed of light. When more 
than two instruments participate in the coincidence all 
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FIG. 5. Examples of joint PDFs for SNRs in the three pairs 
of instruments in the test analysis. 


time differences between pairs of events in the coinci¬ 
dence must satisfy the respective constraints. 

There is a maximum rate at which events can be 
recorded (for example, limited by the sample rate of the 
time series), but this limit is much higher than the mean 
rate and so in the absence of signals the number of events 
recorded in some interval from one filter is adequately 
approximated as a Poisson process. After measuring the 
mean rate of events in a filter in each instrument, say by 
counting the number collected over some large interval 
of time, it is straightforward to use the single-instrument 
rates and the coincidence windows to compute the rates 
at which coincidences will occur between pairs of instru¬ 
ments. For example, if two instruments yield events at 
mean rates fii and ^, 2 , and the coincidence window is ri 2 , 
the mean rate at which coincidences occur is 

M1A2 = (6) 

because fj,i is the number of events in instrument 1 per 
unit time, and /i 2 x 2 ri 2 is the number of events in in¬ 
strument 2 one expects to find within ±ti 2 of each of 
them. 

Computing the rate at which higher-order coincidences 
occur is non-trivial, and we do so using a stone-throwing 
technique. For example, to compute the rate at which 
triple coincidences occur, we start by computing 

M1A2,1A3 = 2^p,i/i2/i3Tl2Tl3) (7) 

giving the rate at which events in instrument 1 are found 
in coincidence with events in both instruments 2 and 3. 
The rate at which 3-way mutual coincidences occur is 

M1A2A3 = f‘iA2,iA3-F’(2 and 3 are comcident|l), (8) 

where P{2 and 3 are coincidentjl) is the probability that 

events in instruments 2 and 3 are found to be in coinci¬ 
dence given that they are both known to be coincident 
with an event in instrument 1. This probability is com¬ 
puted by stone throwing: uniformly-distributed time dif¬ 
ferences are chosen for instruments 1 and 2 and for in¬ 
struments 1 and 3, both consistent with the coincidence 
criteria for those pairs; the time difference between 2 and 
3 is computed from these and if it satisfies the coincidence 
criterion for that pair a successful outcome is recorded; 
the ratio of successful outcomes, n, to total trials, to, 
converges to the desired probability as the number of it¬ 
erations increases. The number of successful outcomes 
is a binomially-distributed RV with standard deviation 
A^/ mp{l — p) < \/ml2 where p is the probability of suc¬ 
cess [3, Section 4.3]. We iterate until the ratio of the 
upper bound of the standard deviation of the number of 
successes to the observed number of successes, n, falls 
below some tolerance 



(9) 


This stopping criterion requires a minimum of e ^/4 it¬ 
erations before it can be met. The number of iterations 
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required, beyond that, is minimized by choosing “instru¬ 
ment 1” to be the instrument for which the Cartesian 
product of the coincidence windows between it and the 
other instruments is smallest, so that coincidence with 
that instrument provides the tightest prior constraint on 
the time differences among the other instruments, and 
the rate of successful outcomes is maximized. 

When some combination of instruments yields mutu¬ 
ally coincident events, it is possible those events are also 
coincident with an event in some other instrument. The 
rate at which coincidences occur that involve some com¬ 
bination of instruments and no others can be found by 
subtracting from that instrument combination’s coinci¬ 
dence rate the rates of coincidences of all proper super¬ 
sets of that combination. For example, if the network 
consists of three instruments, the rate at which exactly 
instruments 1 and 2 are coincident is 


— M1A2 — Ai/1A2A3) (10) 

where ^j,^ia 2 a 3 = /^ia 2 a 3 because there are no other in¬ 
struments. Having obtained, for each instrument com¬ 
bination, the rate at which coincidences occur involving 
those instruments and no others, the probability that 
a coincidence chosen at random involves exactly some 
combination of instruments is obtained by the ratio of 
each such rate to the sum. For example, in the three- 
instrument case, 


M/1A2 


P ({1, 2} jnoise) = - 

+ M/1A3 + M/2A3 + M;!^1A2A3 

( 11 ) 

An example of the result with e = 10 and a compari¬ 
son of the predicted to the observed behaviour is shown 
in FIG. 3. 


E. P (pinst,xLt|noise) 

This PDF is obtained exactly as described in [12], by 
histograming the properties of single-instrument events 
that fail to be coincident with events in any other instru¬ 
ments. As before, the implementation uses a histogram 
binned on a rectangular grid in vs. p in order to 

better ht the natural shapes of the PDF’s contours to the 
rectangular grid of bins. For both the p and x^/p^ sxes 
we now use the tan“^ In binnings described in Appendix 
B. 


F. C(6) = P (P|signal) /P (P]noise) 

In considering the Bayes factor for the intrinsic pa¬ 
rameters of a candidate, 9, we first point out that in 
the example implementation, a “coincidence” is formed 
when two or more instruments in the network register 
events from the same filter within some window of time. 
Therefore, there is no distinction between the 9 carried 


by single-instrument events and the 9 carried by a coin¬ 
cidence. Searches for GWs using algorithms that permit 
a waveform mismatch between instruments in a coinci¬ 
dence will require a different technique for computing the 
Bayes factor for 9 than presented here. 

The denominator, P(djnoise), is obtained as a by¬ 
product of the procedure in Section HID for obtaining 
P({H1,L1,...}|0, noise). Recall that for brevity we have 
been omitting 9 from the parameters of the PDFs, and 
that (11) was implicitly for a choice of 9. The denomina¬ 
tor in (11) gives the total rate of background coincident 
events expected for a choice of 9. P (djnoise) is obtained 
by dividing that rate by the sum of those rates for all 9. 

P (dj signal) is chosen by the individuals performing 
the search. Specifying the relative frequency at which 
templates are expected to yield candidates as a result of 
genuine signals is how the ranking statistic is informed of 
ones prior belief about the distribution of the intrinsic pa¬ 
rameters carried by signals of astrophysical origin. In the 
example implementation demonstrated here, P (dj signal) 
was chosen to be uniform on the manifold in waveform 
space comprising the template bank, i.e., all templates 
are considered equally likely in the signal population. 
This prior is the one implicitly chosen in past searches for 
GWs from CBGs. Because the waveform space’s metric 
is determined by the noise spectrum of the antennas, i.e., 
not by the physical properties of the sources, this prior 
is not uniform in any astrophysically meaningful quan¬ 
tities, like the mass of the source. Other choices might 
be better in future searches. One could conceivably fix 
P (dj signal) to be uniform in template index, and vary 
the template density to affect the desired prior in astro- 
physical parameters. This is the most computationally 
efficient as it avoids, maximally, the filtering of templates 
whose events are subsequently discarded by the ranking 
statistic, but one needs to be certain of the choice of prior 
in advance because the SNR lost in low-density regions 
of the template bank cannot be recovered without re¬ 
analyzing the data if one decides that a different prior is 
desired. Although searches for GBGs have not done so, 
template-based searches for cosmic string bursts [10] and 
for continuous waves (GWs) from pulsars [23] do incor¬ 
porate astrophysical priors into their template weighting 
and/or placement. 


G. Time Dependence 

As discussed above, we have constructed a framework 
that allows us to include the time dependence of the 
antenna noise floors (encoded via horizon distance) in 
the ranking statistic. The daily cycle of human activ¬ 
ity, weather, and thermal changes to optics can alter the 
Gaussian noise floor of the antennas substantially on time 
scales of hours, these changes are easily tracked on time 
scales of minutes, and so there is both a desire and an 
opportunity to incorporate knowledge of Gaussian noise 
floor fluctuations in the assessment of candidates on an 



event-by-event basis. 

Other quantities change in time as well. The mean 
background event rates change, and the antennas un¬ 
dergo periodic maintenance that can alter the statistical 
properties of their noise processes, for example changing 
the distribution of parameters. We do not explicitly 
include the time dependence of these other quantities in 
the description of the ranking statistic given here. In 
practice it is found that a week or more of data is re¬ 
quired before the histograms used to model the noise 
PDFs in the ranking statistic’s denominator have con¬ 
verged enough that the ranking statistic value assigned 
to a given candidate is, effectively, stable. Therefore it is 
not possible to monitor or even define changes to these 
PDFs on time scales less than 0(1) week. 

At the same time, to simplify the management of the 
computer systems, one tends to analyze data in blocks of 
time, and these blocks tend to be of about the same du¬ 
ration as is required to stabilize the denominator PDFs. 
Therefore, through this piecewise analysis of the data, 
by simply collecting fresh histograms within each block 
of time the time dependence of, for example, the p, x^ 
histograms, is naturally incorporated into the ranking 
statistic. Finally, as was done in Section HID to nor¬ 
malize the ranking statistic across different instrument 
combinations, recording the mean background event rate 
in each interval of time allows the ranking statistic to 
be normalized across disjoint analysis blocks providing a 
single, unified, ranking statistic for the entire experiment. 

IV. FALSE-ALARM PROBABILITY 

We now have all the ingredients of the ranking statistic, 
and so at the end of a search for GWs from CBCs we can 
order the candidate events from most signal-like to least 
signal-like. We wish to know how significant a discovery 
the most signal-like of the candidates are: should we have 
expected to see things so signal-like by random chance 
from the noise, or might these be genuine signals? To 
accomplish this, we require 

P(£|0, noise) = f P(... |0, noise) d"“^S, (12) 

the PDF for C (for a choice of intrinsic parameters) 
in the noise population. This is obtained by integrat¬ 
ing the noise PDF for all n parameters in the ranking 
statistic over n — 1 dimensional surfaces, S, of constant 
C. The complementary cummulative distribution func¬ 
tion (CCDF) of (12) gives the probability that a noise 
event drawn at random is found to have a value of C at 
least as high as some threshold, and from that one can ob¬ 
tain the probability that a signal-free experiment of some 
duration (that yields some number of events) yields one 
or more events with values of C at least as high as some 
threshold — the FAP. The details of this were discussed 
in [12], here we will describe how (12) is evaluated in the 
context of the new expression for C. 


We do not know the surfaces of constant £, instead we 
construct an approximation of P(£|0, noise) by visiting 
points in the n-dimensional parameter space, evaluating 
C and P{... |0, noise) at each point that we visit (“...” 
are the n-dimensional co-ordinates of the point), and con¬ 
structing a histogram of C with each sample weighted by 
P(... |0, noise). As explained in Section I, in [12] the 
approximation of the LR as a product of several low¬ 
dimensional terms allowed this sampling procedure to 
be performed by exhaustively visiting every point in the 
discretely-sampled PDFs. Because, here, the numerator 
of C no longer factors into independent single-instrument 
terms we can no longer use this technique, exactly. In¬ 
stead of factoring C and exhaustively evaluating it at each 
grid point in the single-instrument PDFs, we leave it as 
an n-dimensional function and evaluate it at a randomly- 
selected subset of grid points in the n-dimensional space, 
and histogram those results alone. As the number of 
samples gets large the histogram of C samples converges 
to P(£|0,noise). The rate of convergence can be in¬ 
creased by adjusting the PDF from which grid points 
are drawn (and re-weighting the samples appropriately). 
This technique for approximating integrals is known as 
importance-weighted sampling [24, Section 11.7.2]. 

The total number of distinct n-dimensional grid points 
in the binning used to represent C is enormous — in our 
example it’s over 5.6 x 10^^ — but we find that a satisfac¬ 
tory approximation of the PDF can be obtained after just 
4 X 10^ samples. The result is shown in FIG. 6. This was 
obtained using an tan“^ binning for ln/1 as described in 
Appendix B with Inxio = 0, Inxhi = 110, and n = 3, 000. 
All indexes for the binnings comprising the space over 
which C is defined were drawn from uniform distributions 
except the p bins whose indexes were drawn from a distri¬ 
bution whose cummulative distribution function (CDF) 
is proportional to (bin index)0 ®/(# instruments)^* j.g_ 
stricted to bins not below the single-instrument SNR 
threshold, p*. This distribution was arrived at empiri¬ 
cally to efficiently favour interesting SNRs. The counts 
in the \tiC histogram were convolved with a Gaussian 
density estimation kernel having a standard deviation of 
8 bins, the array of counts then rescaled so that its sum 
was 1, and then each bin divided by its size in ln>C to 
yield a normalized, discretely-sampled, PDF. 

Note that the PDF and CCDF have been truncated at 
ln/1 = 5. The procedure by which candidate events are 
collected employs a clustering step to remove redundant 
candidates resulting from each signal. It has been found 
that above In £ = 5 the noise process yields candidates at 
a sufficiently low rate that their population is unaffected 
by the clustering. Modelling the PDF for the noise pro¬ 
cess below ln£ = 5 requires the clustering process to be 
taken into consideration, and this will be the subject of 
future work. The remainder of the procedure for convert¬ 
ing P(£|0, noise) into a mapping from L to FAP and FAR 
is exactly as in [12], and without the clustering model is 
applicable only above the cut-off threshold of ln£ = 5. 
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\nC 

FIG. 6. (Top) PDF for \nC obtained from the importance- 
weighted sampling integration of the noise model. This is, 
in effect, a normalized histogram of the samples constructed 
with the tan“^ In binning described in the text. Note the 
enormous dynamic range on both axes that has been captured 
by the binning and importance-weighted sampling. (Bottom) 
CCDF for the noise model (solid) obtained by integrating 
the PDF in the top panel, shown together with the actual ob¬ 
served CCDF (dashed). The departure of the observed CCDF 
from the model is due to an astrophysically realistic popula¬ 
tion of injections that was present in this data set, but note 
the excellent agreement in the bulk. 


V. EVENT RATE ESTIMATION 

Since we have a reasonable model for P(... |signal) 
in the numerator of £, the procedure described above 
for approximating (12) can also be used to approximate 
P(£|0, signal) by replacing P(... |0, noise) in the inte¬ 
grand with P{.. .\9, signal). In fact, the exact same pro¬ 
cedure is used, and approximations of both P{C\9, noise) 
and P(£|d, signal) are obtained simultaneously in the 
same sampling loop. Following through with a marginal¬ 
ization over 9 yields P(>C|signal). An example is shown 
in FIG. 7, where both the PDF and CCDF have been 
truncated at InP = 5 to match what was done with the 
noise model. 

P(£ I signal) and P(£| noise) are the ingredients re¬ 
quired to implement the signal rate estimation technique 
of Farr et al. [15]. Farr et al. derive the joint PDF for 
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InT 

FIG. 7. (Top) PDF for InT obtained from the importance- 
weighted sampling integration of the signal model. (Bottom) 
CCDF for the signal model obtained by integrating the PDF 
in the top panel. The fluctuations observed in the PDF at 
high values of the ranking statistic appear to be due to sam¬ 
pling noise in the SNR PDFs at high SNR, in particular the 
fluctuations are not the result of under-sampling the rank¬ 
ing statistic: running the stochastic integration for more it¬ 
erations, changing the random number generator’s seed, and 
changing the distribution from which samples are drawn in 
the integral are not found to alter the appearance of the peaks 
and troughs in this PDF, but they are changed by generating 
new SNR PDFs and increasing the number of iterations in 
the SNR PDF sampler decreases their amplitude. 


the rates of signal events and noise events from the rank¬ 
ing statistic values assigned to all events collected in the 
experiment, [15, equation (21)]. Marginalizing this PDF 
over background rate yields the posterior PDF for the 
rate of signals during the experiment. Farr et al. imply 
the possibility of obtaining a closed-form expression for 
this posterior, but for experiments with a large number of 
candidate events we find that approach to be intractable. 
We resort to Markov chain Monte Carlo sampling from 
the joint PDF for the marginalization. For this we use 
the emcee sampler by Foreman-Mackey et al. [25]. 

We wish to obtain credible intervals from the signal 
rate posterior PDF. In particular, we wish to know if 
the 99.9999% credible interval excludes 0. To achieve 
this, we need to measure the posterior’s tails well, and 
running the sampler long enough to do so using the cor- 
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Event Rate Posterior 



FIG. 8. Signal rate posterior PDF. From darkest to lightest 
the shaded areas indicate the 68%, 95% and 99.9999% {“5a”) 
credible intervals. Note that the latter excludes 0. 


rect PDF is inconvenient. Therefore, we sample from 
the square root of the joint PDF given in [15, equation 
(21)], and correct the histogram of samples by squar¬ 
ing the count in each bin and dividing by the bin size. 
An example of the result is shown in FIG. 8. This was 
obtained with 40 walkers, the chain burned-in for 1,000 
iterations, then run for 400,000 iterations, and the sam¬ 
ples histogrammed using a uniform-in-logarithm binning 
spanning the range of values returned from the sampler 
and having about 22,000 bins (so there is an average of 
about 730 samples per bin). Implementing the posterior 
function in C and parallelizing the sum within it using 
OpenMP, the sampling process takes several hours on a 
modern 16-core machine. After correcting the bin counts 
as described above, the corrected counts were convolved 
with a Gaussian density estimation kernel having a stan¬ 
dard deviation of 5 bins, the convolved bin counts nor¬ 
malized to have a sum of 1, and finally divided by the bin 
sizes to yield the normalized discretely-sampled estimate 
of the posterior PDF. 

The mean of the signal rate posterior PDF in the 
example is 8.0 experiment”^, the maximum likelihood 
signal rate is 6.7 experiment”^, the 68% credible inter¬ 
val is [3.6,10] experiment”^ , the 95% credible interval 
[1.8,15] experiment”^, and the 99.9999% credible inter¬ 
val is [.040,38] experiment”^ which, in particular, ex¬ 
cludes 0. In our example analysis, we have the benefit of 
knowing the list of synthetic signals that were added to 
the data, and we can consult that list to see how many of 
them are expected to be visible above the SNR threshold 
in at least two instruments. That number as a function of 
the SNR threshold is shown in FIG. 1. The example anal¬ 
ysis used a single-instrument SNR threshold of p* = 4, 
and apparently 19.4 signals are expected to have yielded 
candidate events involving a coincidence between at least 
two instruments. Additionally, as discussed above, to 
remove the need to model clustering survival, a rank¬ 
ing statistic cut has been imposed discarding all events 
with ln£ < 5. Of the 19.4 signals expected to yield 


coincident events, only 11.5 are expected to survive the 
ranking statistic cut.^ This rate of signals is within the 
75% credible interval, and so there is already reasonable 
agreement between the posterior PDF and the known 
expected signal rate, but it is also possible that a bias 
is present. In considering the origin of a bias, our sim¬ 
ple choice of P(0|signal) = const does not correspond 
to the probabilities with which signals are expected to 
be recovered by different templates, and some bias is to 
be expected as a result of this. Probably the dominant 
reason for a bias, however, is that the injected signals 
were drawn from a population of sources with spinning 
components whereas the templates used to identify can¬ 
didates were non-spinning and therefore would not have 
recovered as much SNR as they could have, and even 
small losses can significantly impact the expected rate. 
For example, assuming the mismatch due to spin results 
in an average loss of SNR of just 6% reduces the expected 
signal recovery rate to 8.9 experiment”^, which is within 
the 50% credible interval, and in essential agreement with 
the posterior PDF. A detailed investigation of the SNR 
recovery efficiency of the template bank used in the ex¬ 
ample search is beyond the scope of this work; the naive 
assumption that the template bank efficiency is given by 
the template bank’s density already yields good agree¬ 
ment between the posterior PDF for the signal rate and 
the signal population, and so if there is a bias in the sig¬ 
nal rate posterior it is not greater than can be explained 
by assuming a few percent loss of SNR. 


VI. CONCLUSION 

We have shown the construction of a new ranking 
statistic for searches for GWs from GBCs, and demon¬ 
strated its implementation including an exposition of the 
technical details required to do so. We have shown that 
the ranking statistic continues to allow us to construct a 
mapping to FAP without the use of time slides. We have 
shows that the ranking statistic can now be used with 
the rate estimation technique of Farr et al. [15] to obtain 
a reasonable posterior PDFs for the rate of signals ob¬ 
served during an experiment. We have found that while 
the identification of single, statistically-significant, out¬ 
lier is not adversely affected by small SNR losses arising 
from signal/template mismatch, the inference of a signal 
rate from the population of candidate events collected by 
the analysis can be sensitive to small SNR losses due to 


^ This is arrived at by drawing SNR and values from distri¬ 
butions constructed for each injection, evaluating the ranking 
statistic at those values, and measuring the survival probabil¬ 
ity for each injection by iterating. That analysis indicates that 
12.1 signals are expected to survive all cuts, but at the ranking 
statistic threshold the clustering is already causing about 5% of 
events to be lost so we suppose that 0.6 more events are lost due 
to clustering and conclude that about 11.5 should remain. 
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template bank inefficiency. This is not unexpected since 
the rate of detectable signals should scale cubically with 
the SNR collected by the templates, and we showed how 
adjustments to the expected rate of recovered signals re¬ 
sulting from as little as a 6% SNR loss in the template 
bank can bring the measured signal rate posterior PDF 
into good agreement with the population of signals we 
know to be in the test data. 

Future work will look more carefully at the problem of 
interpreting the signal rate posterior PDF in terms of an 
astrophysical merger rate. In particular, we will exam¬ 
ine the problem of modelling the loss of low-significance 
events due to clustering in more detail, and incorporate 
a more detailed analysis of the SNR recovery efficiency 
of template bank in the interpretation. 
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Appendix A: Distance and SNR 

For convenience, we collect here some relationships re¬ 
lating the distance to a GW source, its SNR, and the 
number density of sources. 

The strain seen in a GW antenna is the projection of 
the GW field onto the antenna’s response tensor, and in 
terms of the two transverse, traceless, polarization com¬ 
ponents of the field, h+ and fix, can be written [27, Sec¬ 
tion 9.4] 

h{t) = F+{(j),0,'tp)h+(t) + Fx{(j),e.2p)hx{t). (Al) 

where the antenna response factors, F+ and F ^, depend 
on the direction to the source, {(j), 9), and the orientation 
of the polarization axes, ip [28, equation (B6)]. For non- 
precessing CBGs whose radiation is dominated by the 
l,m = 2,2 mode, h{t) is inversely proportional to the 
effective distance [29, equation (3.3c)] 

_ 1 

DeS = D FF^cos^i , (A2) 


where D is the physical distance to the source, and l the 
angle between the line-of-sight to the source and its or¬ 
bital axis. Deff > D, and a source is “optimally oriented” 
with respect to an antenna when Dgg = D. 

Defining po to be the nominal matched-filter output 
observed for a signal in the absence of noise, po oc h{t) oc 
D~g [29]. This relationship can be written as 

Po = SDn/Deff, (A3) 

thereby defining Dh, the “horizon distance”, the physi¬ 
cal distance at which an optimally-oriented source is seen 
with a nominal SNR of 8 [29]. Note that Dh depends on 
the antenna noise spectrum and the physics of the GW 
source: one may speak of an horizon distance being as¬ 
sociated with a source by assuming a canonical antenna 
noise spectrum, or with an antenna by assuming a canon¬ 
ical source; stronger emitters of GWs are said to have 
larger Du, less sensitive antennas smaller Du. Combin¬ 
ing (A2) and (A3), 

1 

Dpo = 8Dh FF^cos^i = D, 

(A4) 

where the parameter D is introduced for compactness: 
a sort-of “noise-limit distance”, the physical distance at 
which a source with the given geometric arrangement 
with respect to a given antenna is seen at an SNR of 
1 in that antenna. 

In stationary Gaussian noise the SNR, p, at which a 
signal is recovered is a 2 degree-of-freedom (DOF) non¬ 
central x-distributed RV with noncentrality parameter 
Po [29, Section IV]. This is also known as the Rice distri¬ 
bution with (7 = 1 and noncentrality parameter po [30, 
equation (3.10-11)].^ The probability that a source with 
nominal SNR po is observed above some SNR threshold 
p* is 

P{p > P*\po) = Qi{po,P*) = Qi(^,p*), (A5) 

where Qi is the first-order Marcum Q function [31]. As¬ 
suming the number density of sources to be uniform in 

volume, it is oc D"^ AD = —P^Apo, therefore the total 

^ 0 . . . 

number of sources in a given direction with a given po¬ 
larization and orbit inclination visible to an antenna with 
a given horizon distance is 

nOO 

# sources (xD^ Qi (/?o, dpo. (A6) 

Jo 

Because Qi{0,p*) ^ 0 the integral diverges. Physically, 
assuming a uniform source density to arbitrary distance 


^ Reference is to Rice’s original work; he does not call the distri- 
bution by that name. 
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yields an infinite number of sources at infinite distance 
with zero nominal SNR that, nevertheless, due to noise 
fluctuations have a non-zero probability of registering as 
events. Because the universe has a finite age and stellar 
evolution must be allowed to progress for a period of 
time after the birth of the universe before compact object 
mergers can occur, the source density must, really, fall to 
0 at some distance, and if this is correctly accounted for 
this integral must be finite — Olbers’ paradox revisited. 

If the number density cutoff occurs at a distance D ':$> 
D (the antennas cannot see as far as the first sources), 
then the integral in (A6) becomes a function only of the 
signal detection threshold, p*, and taking that to be a 
fixed parameter we can side-step the divergence and say 


every case the variable is valid from 0 to oo, and so an in¬ 
finite number of uniform-in-the-logarithm bins would be 
required to define the function everywhere. To address 
this, we make use of bins that are uniform in the arct¬ 
angent of the logarithm of the variable. The arctangent 
function is approximately linear near 0 and so by choos¬ 
ing an appropriate translation and scaling one obtains a 
binning that is approximately logarithmic over some cho¬ 
sen range of values, but that spans a domain from 0 to 
-boo with a finite number of bins. 

An tan“^ In binning can be defined by three parame¬ 
ters: x\o and Xhi, the low and high roll-offs defining the 
range of approximately logarithmic binning, and n, the 
total number of bins. The bin boundaries are given by 


# sources oc D^. (A7) 


Appendix B: tan ^ In Binning 


Xk = exp 


'c-2 

/"K , 7r\ 

d — tan 

n) 

TT 

Kn 2 / 


-b Ina; 


0 < k < n, 

(Bl) 


In implementing the ranking statistic described in this 
article, one will need to histogram randomly-drawn sam¬ 
ples that span a large domain and whose density varies 
greatly over that domain. It is helpful to use non-uniform 
binnings whose bin density is approximately inversely 
proportional to the sample density so that the number 
of samples falling in each bin — and thus the error from 
counting fluctuations — is approximately constant. 

The functions described by the binned samples that 
are encountered here tend to be well described by sim¬ 
ple polynomials in the logarithm of the function and the 
logarithm of the variable (see, for example, any of the 
functions depicted in FIG. 5, FIG. 6, or FIG. 7), and so 
binnings that are uniform in the logarithm of the variable 
would seem to be convenient. Unfortunately, in almost 


where lux = (lna;hi+lnxio)/2, and 6 = (Inxhi —lna:io)/2. 

One must be careful evaluating this function numeri¬ 
cally, but, still, except for binnings with very few bins, 
one nearly always finds that underflows and overflows 
prevent the representation of the first few and last few bin 
boundaries using double-precision floating-point num¬ 
bers, making the first few and last few bins appear to 
have identical upper and lower boundaries. We simply 
discard those bins, therefore generally the number of use¬ 
ful bins is slightly less than n. 

For example, if a;io = 10, xm = 100, and n = 10, the 
bin boundaries are 0, 3.314, 11.53, 18.57, 24.92, 31.62, 
40.13, 53.86, 86.72, 301.8, oo. The factors separating 
adjacent boundaries are oo, 3.48, 1.61, 1.34, 1.27, 1.27, 
1.34, 1.61, 3.48, oo. 
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